Hyperbolic slicings of spacetime: singularity avoidance and gauge shocks 
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I study the Bona-Masso family of hyperbolic slicing conditions, considering in particular its prop- 
erties when approaching two different types of singularities: focusing singularities and gauge shocks. 
For focusing singularities, I extend the original analysis of Bona et. al. and show that both marginal 
and strong singularity avoidance can be obtained for certain types of behavior of the slicing con- 
dition as the lapse approaches zero. For the case of gauge shocks, I re-derive a condition found 
previously that eliminates them. Unfortunately, such a condition limits considerably the type of 
slicings allowed. However, useful slicing conditions can still be found if one asks for this condition 
to be satisfied only approximately. Such less restrictive conditions include a particular member of 
the 1+log family, which in the past has been found empirically to be extremely robust for both Brill 
wave and black hole simulations. 



PACS numbers: 04.20.Ex, 04.25. Dm, 95.30.Sf, 



I. INTRODUCTION 

The choice of good coordinates plays a crucial role 
in finding solutions of the Einstein equations. This is 
particularly important in the case of numerical simula- 
tions of strongly gravitating systems, where a bad coor- 
dinate choice can easily lead to the formation of a coor- 
dinate singularity which stops the numerical simulation 
and severely limits the region of the spacetime covered 
by it. Coordinate singularities are not the only concern: 
the presence of physical singularities, such as those asso- 
ciated with black holes, can also have a deep impact in a 
numerical simulation, as unless special care is taken the 
time slices can march right onto the physical singularity 
very early during a simulation. 

When considering dynamical evolutions of spacetime 
based on a 3+1 decomposition [1, 2], the coordinate 
choice naturally separates in two different aspects: the 
choice of a specific foliation of spacetime into spatial 
hyper-surfaces (also known as the "slicing"), associated 
with the lapse function a, and the choice of the way in 
which the lines of constant spatial coordinates (the "time 
lines") propagate from one hyper-surface to the next, as- 
sociated with the shift vector (3 % . Here I will concentrate 
fully in the role played by the choice of a slicing condition 
and leave the choice of a shift vector for a later work. 

In order to specify a foliation of spacetime, one needs to 
prescribe a way to calculate the lapse function a, which 
measures the proper time interval between neighboring 
hyper-surfaces along their normal direction. There is, 
of course, an infinite number of ways in which one can 
choose the lapse function, but typically the different 
choices can be classified in the following way: 1) Pre- 
scribed slicings, where the lapse is specified as an a pri- 
ori known function of space and time, 2) algebraic slic- 
ing conditions, where the lapse is specified algebraically 
as some function of the geometric variables (metric and 
extrinsic curvature) at each hyper-surface, 3) elliptic slic- 
ing conditions, where the lapse is obtained by solving an 
elliptic differential equation at every time step that typ- 



ically enforces some geometric condition on the spatial 
hyper-surfaces, and 4) time derivative slicing conditions, 
where the time derivative of the lapse is specified as some 
algebraic function of the geometric quantities and the 
lapse is evolved as just another dynamical quantity (this 
last case is often included in the algebraic slicing class 
mentioned above). 

An example of a prescribed slicing is the so-called 
"geodesic slicing" , where one simply takes a — I. An 
example of an algebraic slicing condition is "harmonic 
slicing" where one takes instead a — */7, with 7 the de- 
terminant of the spatial metric. A well known elliptic 
slicing condition is the "maximal slicing" condition [3], 
which requires that the spatial volume elements remain 
constant during the evolution. Elliptic conditions are 
typically robust and well behaved, but have the drawback 
of being computationally expensive. Algebraic conditions 
are much easier to apply, but are difficult to analyze in a 
general case. Time derivative slicing conditions, on the 
other hand, have the advantage of being both easy to 
implement and, in the particular case when they lead to 
hyperbolic equations (as is the case of the Bona-Masso 
family [4], see below), much easier to analyze as well. 
They also include many well known algebraic conditions 
as their integral form. 

In this paper, I will only consider hyperbolic slicing 
conditions and study in which ways such conditions can 
lead to pathological slicings and how can those patholo- 
gies best be avoided. There are many different ways 
in which a foliation of spacetime can become patholog- 
ical: The slices can hit a physical singularity, the slices 
can hit a coordinate singularity where the volume ele- 
ments become zero (the normal lines focus), the slices 
can become non-smooth at some point or the slices can 
remain smooth but stop being space-like (they can be- 
come null at a point, for example). Of the different pos- 
sible pathologies mentioned above, I will concentrate on 
two specific types: "focusing singularities" [5] defined as 
those for which the spatial volume elements vanish at 
a bounded rate, and "gauge shocks" [6] defined as so- 
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lutions for which the lapse becomes discontinuous as a 
consequence of the crossing of the characteristic lines as- 
sociated with the propagation of the gauge, and the time 
slices therefore become non-smooth. 

This paper is organized as follows: In Sec. II, I intro- 
duce the Bona-Masso hyperbolic slicing condition. Fo- 
cusing singularities are defined in Section III where I also 
find under which circumstances the Bona-Masso slicing 
condition is singularity avoiding. In Sec. IV, I introduce 
the idea of a gauge shock, derive a condition to avoid 
them and see what slicings obey that condition cither 
exactly or approximately. I conclude in Sec. V. 



II. THE BONA-MASSO FAMILY OF 
HYPERBOLIC SLICING CONDITIONS 

The Bona-Masso family of slicing conditions [4] is a 
time evolution type of condition for which the lapse is 
chosen to satisfy the following evolution equation 
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a=(d t -C )a = -a 2 f(a)K, 



(1) 



with Cp the Lie derivative with respect to the shift vec- 
tor P l , K the trace of the extrinsic curvature and f(a) a 
positive but otherwise arbitrary function of a. The rea- 
son why f(a) has to be positive will become clear below. 
Here we just mention the fact that the right hand side of 
condition (1) is the most general term one can construct 
that involves only first order spatial scalars [4]. 

There are several things that are important to men- 
tion about the family of slicing conditions (1). First, 
we notice that even if this family was proposed in the 
context of the Bona-Masso hyperbolic reformulation of 
the Einstein equations [4, 7, 8, 9, 10], it is in fact quite 
general and can be used successfully with any particular 
form of the 3+1 evolution equations. This was shown 
recently when such a condition was used together with 
the Baumgarte-Shapiro Shibata-Nakamura (BSSN) for- 
mulation [11, 12] to obtain long-term stable and accu- 
rate evolutions of black hole spacetimes [13, 14]. Also, 
condition (1) is a generalization of slicing conditions that 
have been used in evolution codes based on the Arnowitt- 
Deser-Misncr (ADM) formulation [1, 2] since the early 
90's [15, 16]. 

Second, condition (1) can also be trivially adapted to 
the case when, instead of the lapse a, one evolves a den- 
sitized lapse of the form 



Q 



a 7 



r/2 



(2) 



with 7 the determinant of the spatial metric 7^ and a 
a constant parameter. Such a densitized lapse (particu- 
larly the case a = — 1) has recently been advocated in 
the context of hyperbolic reformulations of the Einstein 
equations (see for example [17, 18]). In terms of Q, con- 
dition (1) becomes 



(3) 



where to calculate the Lie derivative with respect to j3 % 
contained in the operator d/dt one must use the fact that 
Q is a density of weight a. 

Finally, it is also important to mention that the shift 
terms included through the Lie derivative in condition (1) 
are such that one is guaranteed to obtain precisely the 
same spacetime foliation regardless of the value of the 
shift vector. This would seem to be a natural require- 
ment for any slicing condition. However, it is plausible 
that in a particular situation one would like to choose a 
slicing condition and a shift vector that are closely inter- 
related (see for example [19]). Indeed, generalizations 
of the Bona-Masso slicing condition that in the pres- 
ence of a non-zero shift vector do not have the form (1) 
have already been used in the literature. For example, 
Refs. [20, 21] use d t a = -af(aK - D i f3 i ) instead of (1) 
(with Di the spatial covariant derivative), and Ref. [14] 
simply uses d t a = -a 2 f(K - K ) (with K = K(t = 0)) 
even with a non-zero shift vector. Which off these gen- 
eralizations is best under different circumstances is an 
important question that I will leave for a future work. 

With these comments in mind, let us go back to con- 
dition (1). Taking an extra time derivative we find 



d 2 2 , 
d^ a = - af 



-K a(2f + af)K 2 



(4) 



with /' := df /da. From the ADM evolution equations 
one easily finds that, in vacuum, 



^-K = -D 2 a + aK ij K ij , 
dt 



(5) 



with D 2 the Laplace operator associated with the spatial 
metric. This last equation implies 



dt 2 



a — a fD a = 

~a 3 f[K lJ K^-(2f + a.f')K 2 ] 



(6) 



Equation (6) shows that the lapse obeys a wave equa- 
tion with a quadratic source term in . It is because of 
this that we say the slicing condition (1) is a hyperbolic 
slicing condition: it implies that the lapse evolves with 
a hyperbolic equation (but see the discussion on gauge 
shocks below for a more formal proof of hyperbolicity). 
The wave speed associated with equation (6) along a spe- 
cific direction x % can be easily seen to be 



(7) 



Notice that this will only be real if f(a) > 0, which 
explains why we asked for f(a) to be positive. In fact, 
/(a) must be strictly positive because if it was zero there 
would be no complete set of eigenvectors and we would 
not have a strongly hyperbolic system (see Sec. IV B). 

To see how the gauge speed v g is related to the speed 
of light consider for a moment a null world-line. It is 
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not difficult to find that such a world-line will have a 
coordinate speed along direction x % given by 



(8) 



so the gauge speed (7) can be smaller or larger that the 
speed of light depending on the value of /. 

Notice that having a gauge speed that is larger than 
the speed of light does not lead to any causality viola- 
tions, as the superluminal speed is only related with the 
propagation of the coordinate system. One could argue 
that superluminal gauge speeds are not desirable as they 
would allow gauge effects to propagate out of black hole 
horizons, for example. Empirically, however, the most 
successful slicing conditions for the simulation of black 
hole spacetimes have been precisely those that have su- 
perluminal (even extremely large) gauge speeds: maxi- 
mal slicing, which as an elliptic condition has, at least 
formally, an infinite speed of propagation, and the 1+log 
slicing condition which is a member of the Bona-Masso 
family that has superluminal gauge speeds whenever the 
lapse is small (see following section). 



A. Relating lapse and spatial volume elements 

The ADM evolution equation for the spatial metric 7^ 
is given by 



d „ Tjr 

— 7ij = -2aKij , 



(9) 



which implies the following evolution equation for the 
spatial volume elements 7 1 / 2 : 



dt 



r 



/2 



-cry 



l l 2 K. 



(10) 



Comparing this with equation (1), and taking / = 1, 
we can trivially solve for a in terms of 7 1 / 2 to obtain (in 
the case of zero shift vector) 



(11) 



with h{x l ) a time independent function, ft is very im- 
portant to stress the fact that the previous relation holds 
only when moving along the normal direction to the hy- 
persurfaces, and not when moving along the time lines 
which will differ from the normal direction for any non- 
zero shift vector. That is, we are relating the lapse to 
the volume elements as seen by the normal observers. In 
the following, whenever we relate the lapse to the vol- 
ume elements, it should always be understood that we 
are referring to the volume elements associated with the 
normal observers. 

It is not difficult to show (see following section) that 
equation (11) is equivalent to the condition 



™ = <r C = 



(12) 



with g^y the spacetime metric. That is, / = 1 corre- 
sponds to the case when the time coordinate is a har- 
monic function. Because of this the case / = 1 is known 
as "harmonic slicing". Notice also that in this case the 
gauge speed is identical to the physical speed of light, i.e. 
the gauge propagates along null lines. Harmonic slicing 
can be seen to be equivalent to having a time indepen- 
dent lapse density Q of weight a = — 1 (again, assuming 
the shift vanishes). 

One can construct other well known families of slic- 
ing conditions by choosing different forms of f(a). For 
example, if we choose f(a) — N, with N a constant, 
we obtain what can be called the "generalized harmonic 
slicing condition" , which can also be easily integrated to 
give 



a = h(x l ) 7 



N/2 



(13) 



And if we take f(a) — N/a we obtain the "1+log" 
family [22, 23], which again can be integrated to find 



a = h{x l )+\n{j N ' 2 ^j . 



(14) 



The 1+log family mimics maximal slicing and has 
strong singularity avoiding properties (see section III be- 
low). In particular, 1+log slicing with N = 2 has been 
found empirically to be very robust when evolving black 
hole spacetimes [13, 14, 24]. As mentioned before, one 
can easily see that the gauge speed associated with the 
1+log family can become far larger than th e spe ed of 
light as the lapse becomes smaller (v g /vi — y/N/a). 

More generally, using equation (10) one can find that 
for arbitrary f(a) the following relation between a and 
7 1 / 2 holds 



da 
af(a) 



which implies 



7 



1/2 



F(x l ) exp 



f f da 
[J o~f(a~ 



(15) 



(16) 



with F{x l ) again a time independent function. This last 
expression will be the starting point when we discuss fo- 
cusing singularities below. 



B. The foliation equation 

Consider now a spacetime with metric g^ v , and assume 
that we have a foliation of this spacetime into spatial hy- 
persurfaces. Such a foliation allows us to define a time 
function T whose level sets correspond to the members 
of the foliation. One can show that the Bona-Masso slic- 
ing condition (1) can be written as a generalized wave 
equation for the time function T in the following way 



1 - 



1 



/(«) 



T;„ v = , 



(17) 
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with n M the unit normal vector to the spatial hypersur- 
faces. I will call equation (17) the "foliation equation". 
Notice that for harmonic slicing the above equation re- 
duces to the simple wave equation (12). 

The fact that the foliation equation above is equivalent 
to the Bona-Masso slicing condition (1) can be shown by 
choosing 3+1 coordinates {t,x 1 } adapted to the folia- 
tion. In this coordinate system we can take T = t, so 
equation (17) becomes 



1 



1 



/(«) 



. 



(18) 



Using now the 3+1 expressions for the components of the 
4-metric g^ v we obtain the following expressions for the 
Christoffcl symbols 



1 Oi 



r° 



- (dta + Fdia-FPKy) 
a 

- (d ia - (3 m K im ) , 
a 

--K iA . 



(19) 
(20) 
(21) 



Substituting this into (18), and using the 3+1 expressions 
for the normal vector we find 



d t a - Pd t a + a 2 f(a)K = 



(22) 



which is precisely the Bona-Masso condition (1). 

Notice that if we take / > 1, one can always have a unit 
normal vector n M such that the coefficient of the <9 2 T term 
in the foliation equation changes sign, and the equation 
apparently becomes elliptic. The system is in fact still 
hyperbolic (see Sec. IV B), and the change in signature 
just reflects the fact that for / > 1 the characteristic 
cones can tilt beyond the time axis. 

The foliation equation (17) will prove to be very im- 
portant when we study the formation of gauge shocks. 



III. FOCUSING SINGULARITIES 

A very important property of slicing conditions is that 
of "singularity avoidance" . Singularity avoidance refers 
to the property of certain slicing conditions of slowing 
down coordinate time, by making the lapse go to zero, 
when the spatial volume elements go to zero (this is 
known as the "collapse of the lapse"). Recent advances 
in black hole excision techniques [25, 26, 27, 28] would 
seem to minimize the need for singularity avoidance in 
the choice of slicing conditions. One should remember, 
however, that singularity avoidance is not only important 
when one is interested in studying black hole spacetimes 
where real physical singularities are present, but is also 
needed in order to prevent the formation of coordinate 
singularities caused by the focusing of the normal lines 
in regions with strong gravitational fields. 

Bona et. at have shown [10] that the slicing condi- 
tion (1) can avoid so-called "focusing singularities" for 



some choices of the function /(a). Here I will extend 
their analysis and show explicitly what type of behavior 
f(a) must have as a approaches zero in order to avoid 
such singularities. 

Following [10], we define a focusing singularity as a 
place where the spatial volume elements 7 1 / 2 vanish at 
a bounded rate. Let us assume that such a singularity 
occurs after a finite proper time t s away from our initial 
time slice (as measured by the normal observers). From 
the definition of the lapse we see that the elapsed coor- 
dinate time will then be 



At 



Jo a 



(23) 



We will further characterize the singularity by the rate 
at which 7 1 / 2 approaches zero as a function of proper 
time. We will say a singularity is of order m if 7 1 / 2 
approaches zero as 



7 



1/2 



1 \ r 



(24) 



with m some constant power. Notice that m must be 
strictly positive for there to be a singularity at all, and 
it must be larger than or equal to 1 for the singularity to 
be approached at a bounded rate. 

As the volume elements 7 1 / 2 approach zero, there are 
clearly three possible behaviors for the lapse: 1) a re- 
mains finite as 7 1 / 2 vanishes, 2) a vanishes as 7 1 / 2 van- 
ishes, and 3) a vanishes before 7 1 ' 2 vanishes. Case 1 
would clearly imply that coordinate time remains finite at 
the singularity, so the singularity would not be avoided. 
However, from equation (16) one can easily see that if the 
lapse remains always finite it is impossible for the volume 
elements to ever vanish (remember that f(a) is never 
allowed to be zero). We then conclude that case 1 can 
never happen, which implies that the Bona-Masso slicing 
condition (1) always causes the lapse to collapse when 
the volume elements approach zero, for any f(a) > 0. 
Case 3, on the other hand, implies that the time slices 
stop advancing a finite coordinate time before the singu- 
larity is reached (the time slices can in fact move back 
under certain conditions, see below). We will call such 
behavior "strong singularity avoidance" . Finally, case 2 
corresponds to the case when the lapse becomes zero at 
the same time as the volume elements. Whether in such 
a case the singularity is reached after a finite or infinite 
coordinate time will depend on the speed at which a ap- 
proaches zero at the singularity. We will say that a slicing 
is "marginally singularity avoiding" if the singularity is 
reached after an infinite coordinate time. 

To study under which conditions we can have strong 
or marginal singularity avoidance we must say something 
about the form of the function f(a). From now on we will 
therefore assume that, as a approaches zero, the function 
f(a) behaves as 



f(a) = Aa n , 



(25) 
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with both A and n constants and A > 0. Such an as- 
sumption implies that 



1 f da 



da 

af(a) ~ A J a n+1 
lna^ A 



n = 
-\/{nAa n ) n ^ 



(26) 



Let us first consider the case n ^ 0. From equation (16) 
we now find that 



j 1 / 2 exp 



1 



nAa n 



(27) 



As a approaches zero we have two separate cases depend- 
ing on the sign of n: 

.. wo f finite n < /noX 
hm Y' = { n ^_ n ( 28 ) 

Since for n < the volume elements remain finite as the 
lapse approaches zero we conclude that such a case cor- 
responds to strong singularity avoidance. On the other 
hand, for n > both the lapse and the volume elements 
go to zero at the same time so we can at most have 
marginal singularity avoidance. 

For the case n = we find, again using (16), that 



7V2 ~ a iM 



a ~ 7 



A/2 



(29) 



It is then clear that in this case a and 7 1 / 2 also vanish 
at the same time. 

We now need to decide if the cases n > and n = 0, 
for which the lapse and the volume elements become zero 
at the same time, reach the singularity in an infinite or 
a finite coordinate time. For this we need to study the 
behavior of a as a function of proper time r as we ap- 
proach the singularity. Starting from equation (23) for 
the elapsed coordinate time we find 



At 



Jo at 

-f 



dr/da 



da 



(30) 



where ao is the initial lapse and where we are already as- 
suming that we are interested in the case when a vanishes 
at t s . Equation (30) implies that if dr/da remains differ- 
ent from zero as we approach the singularity then At will 
diverge and we will have marginal singularity avoidance. 
On the other hand, if dr/da vanishes at the singularity 
as a p (with p some positive power) or faster, then the 
integral will converge and the singularity will be reached 
in a finite coordinate time. 

To find the behavior of dr/da as we approach the sin- 
gularity, we notice that equation (24) implies 



dj 1 / 2 /dr ~ — m (r s — t) t 



(31) 



and 



rfln 7 V2 



dr (r s - r) 

From this, together with equation (15), we now find 

da/dr m 
af(a) (r s - r) 

which can be integrated to give 



(32) 



(33) 



r, - exp I — 

m 



da 
af(a) 



(34) 



If we now take f(a) given by (25) we finally find 

/ 1 r da 

\ mA J ma nJt 



{ 



exp 



_ a l/mA 



n = 

r s -exp[-l/(nmAa")] n>0 



(35) 



The case n < is not of interest here since we already 
showed that for such a case the lapse will vanish before 
we reach the singularity. 

Let us consider the case n > first. The derivative of 
r with respect to a then turns out to be 



dr 
da 



1 



mAa n+1 



exp 



1 



nmAa n 



(36) 



from which it is easy to see that as a approaches zero, 
dr/da also approaches zero faster than any power. As 
we have seen, this means that the singularity is reached 
in a finite coordinate time, so the case n > does not 
avoid the singularity (not even marginally). 
For n = 0we have, on the other hand 



dr 
da 

We then see that 



lim — 

a^o da 



1 

mA 



a 



1/mA-l 





-1 

-00 



mA < 1 
mA = 1 
mA > 1 



(37) 



(38) 



The case mA < 1 therefore reaches the singularity in a 
finite coordinate time, while the cases mA > 1 reach it in 
an infinite coordinate time and are therefore marginally 
singularity avoiding. 

Our final result can be summarized as follows: If /(a) 
behaves as / = Aa n for small a and we have a singularity 
of order to, then 

1. For n < we have strong singularity avoidance. 

2. For n — and mA > 1 we have marginal singularity 
avoidance. 

3. For both n > 0, and n — with mA < 1, we do not 
have singularity avoidance, even though the lapse 
collapses to zero at the singularity. 
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In the particular case when we have a singularity of 
order m = 1, then harmonic slicing (n = 0, A = 1) 
marks the boundary between avoiding and reaching the 
singularity. 

As a final observation, let us consider again the case 
n < 0, for which we have shown that we have strong sin- 
gularity avoidance. Looking at our original slicing con- 
dition (1) we see that if n < —2 then, as the lapse ap- 
proaches zero, one can not guarantee that d t a will also 
approach zero. The lapse can therefore easily become 
negative and the slices will not only avoid the singularity 
but can in fact back away from it. This type of behavior 
is probably not desirable, as one runs the risk of having 
the time slices stop being space-like (they advance in one 
region and move back in another). If we want to guar- 
antee that we have strong singularity avoidance without 
the lapse becoming negative we must limit ourselves to 
the region — 2 < n < 0. Notice that the 1+log family 
corresponds to n = —1, and is precisely in the middle 
of this range, which probably accounts for the fact that 
empirically it has been found to be a very good choice. 



IV. GAUGE SHOCKS 

In physics, one talks about "shock waves" as solutions 
to the hydrodynamic equations where very sharp density 
gradients propagate through a medium at speeds that are 
higher than the speed of sound in that medium. Math- 
ematically, shocks are discontinuous solutions of a non- 
linear hyperbolic system of equations characterized by 
the fact that characteristic lines converge at the disconti- 
nuity The discontinuity propagates with a speed called 
the "shock speed" that is somewhere in between the val- 
ues of the characteristic speeds in the regions behind and 
in front of the shock. Usually, in order to completely 
determine the form and speed of a shock one needs to 
supplement the evolution equations with extra conditions 
coming from physical considerations known as "entropy- 
conditions" (see for example [29] ) . Such entropy condi- 
tions are necessary because once a discontinuity forms, 
the possible mathematical extensions through it are no 
longer unique and one needs a mechanism to choose the 
physically allowed solutions. 

It is well known that physical shocks (i.e. shocks in 
the geometry) do not appear in general relativity. So- 
lutions of the Einstein equations normally called "shock 
fronts" refer to discontinuities in the curvature of space- 
time present in the initial data that propagate with the 
speed of light. These type of solutions are not shocks in a 
mathematical sense but are called instead "contact dis- 
continuities" . Here we will not consider discontinuities 
in the geometry but rather solutions to our hyperbolic 
slicing conditions that start from smooth initial data and 
develop discontinuities later, when the characteristic lines 
associated with the gauge cross. Since this is the defining 
property of a shock, we will call those solutions "gauge 
shocks". However, we must stress the fact that, in con- 



trast to the case of hydrodynamics, once a gauge shock 
forms one should make no attempt to continue the so- 
lution any further since such a shock will indicate that 
our coordinate system has broken down, and there is no 
physical principle that can be used to extend the solu- 
tion beyond this point. We will therefore not have shock 
waves as such (i.e. propagating shocks), but just shock 
formation. 

Gauge shocks were first studied in References [6, 30], 
where it was found that discontinuities in the lapse can 
easily develop starting from smooth initial data in a wide 
variety of cases. It was also shown how, in some partic- 
ular cases, one can even predict the exact time when a 
gauge shock would form by just analyzing the initial data. 

In [6] a particular condition on the function f(a) was 
found for which gauge shocks would not form: 

1 - / - af/2 = . (39) 

This condition, however, was derived using a hyperbolic 
formulation of the 3+1 evolution equations (the Bona- 
Masso formulation). In the following sections I will re- 
derive the condition making no reference to the Einstein 
equations in any form. 

A. Linear degeneracy and shocks 

Consider a system of equations of the form 

d t u i + d x F i = q i ie{l,...,JV„}, (40) 

where Fi and qt are arbitrary, possibly non-linear, func- 
tions of the u's but not their derivatives. Notice that the 
system above can be written also as 

d tUl + J2 Mijd xUj = qi i e {1, N u } , (41) 
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with My = dFi/duj the Jacobian matrix. 

Let Aj be the eigenvalues of the Jacobian matrix M. 
The system of equations is called "hyperbolic" if all the 
Xi are real. Further, the system is said to be "strongly 
hyperbolic" if the matrix M has a complete set of eigen- 
vectors ej . Let us assume that this is the case, we then 
define the eigenfields Wi in the following way 

u = Rw => w = iT 1 u, (42) 

where R is the matrix of column eigenvectors e^. One 
can show that the matrix R is such that 

RMR- 1 = A , (43) 

with A = diag(Ai). The evolution equation for the eigen- 
fields Wi then turns out to be 

d t Wi + A 4 d x u>i = q\ , (44) 

with q[ a function of the w's but not their derivatives. 
In terms of the eigenfields the system transforms into a 
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series of coupled advection equations with characteristic 
speeds given by the eigenvalues A,. 

If a given eigenvalue Xi is independent of its corre- 
sponding cigenfield Wi, then we say that the eigenfield is 
"linearly degenerate" [29, 31], that is 



dwi 



Ed Xi duj ^_ . 



duj dw. 
j= i i 



(45) 



Linear degeneracy is a sufficient condition for there 
not to be shocks associated with a given eigenfield. One 
can understand this intuitively by noticing that linear 
degeneracy implies that the characteristic lines do not 
change in response to changes in the field propagating 
along them. 



B. Avoiding gauge shocks 

In order to study the effects of our slicing condition 
without having to worry about the evolution of the space- 
time itself we will now assume that we have a known 
background spacetime with metric g^. In that space- 
time, we will consider some initial spatial slice, and then 
construct a foliation according to our slicing condition. 
Let T be the time function associated with our folia- 
tion (that is, each spatial hypersurface will correspond 
to T = constant). In Sec. II B above we saw that for the 
Bona-Masso family of slicing conditions, the function T 
will satisfy the following foliation equation 



l 



l 



V At V,T = 0, 



(46) 



where is the unit normal vector to the hypersurfaccs 
and where V M denotes covariant differentiation with re- 
spect to the 4-metric <7 M „. The unit normal vector 
can be easily constructed from the time function T in 
the following way 



n,. 



-V„T VT) 



1/2 ' 



(47) 



where the overall minus sign is there to guarantee that 
we have a future pointing normal vector. 

Let us now calculate the increment in T if we move a 
proper distance dr along the normal direction 

dT = {dr n") V P T = dr {-V„T V^T) 1/2 . (48) 

On the other hand, from the definition of a we have 
dr = adT. Comparing both expressions we find 



a=(-V (1 TVTr 1/2 . 
The normal vector then takes the form: 



(49) 



(50) 



Consider now a particular point V in spacetime. To 
study the evolution of T close to that point we construct 
locally flat coordinates (t,x l ), so the metric close to V 
becomes the flat metric r}^ and the Christoffel symbols 
vanish. Equation (46) then reduces to 



']d^d v T = Q, 



(51) 



where we have defined a := l/f — 1. Expanding this 
equation we find 

- (1 + a(n ) 2 ) d 2 T + - onV) d l d ] T 

-2an"n l d t d t T = . (52) 

Let us now define II := d t T and := d{T. Equa- 
tion (52) can then be transformed into the system 



\ - 2an ( V 



+ a{n°) 2 

an l n 
IK)" 2 



ES 13 — an l n 3 
1 + a(n°) 2 dl ' 



dill. 



(53) 
(54) 



In our locally flat coordinate system, the contravariant 
components of the unit normal vector become 



n° = +all , n l = - 
and the lapse reduces to 

a = (n 2 -* 2 )- 1/2 



(55) 



(56) 



with * 2 = £\ *? • The system (53)- (54) then takes the 
following form 



' 1 + aot z ii z 

i 

+ 4^ l + aa 2 n 2 1 3 



an. 



(57) 
(58) 



In order to see if our system of equations is hyperbolic 
we consider derivatives along a fixed spatial direction, say 
x, and neglect derivatives along different directions. It is 
evident that in this case the variables with q ^ x, 
can be considered as fixed and we only need to analyze 
the sub-system {II, ^ x }- The Jacobian matrix M for this 
reduced system becomes: 



2aa 2 U^ x 1 - aa 2 ^ 2 



M 



1 + aa 2 n 2 1 + aa 2 II 2 
1 



(59) 



The eigenvalues of this matrix are easily found to be 
1 



{aa 2 U^ a 



1 + aa 2 n 2 

±[i + aa 2 (n 2 -* 2 )] 1/2 } , (60) 
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with corresponding eigenvectors 

e± = (A±,l) . 



(61) 



Having found the eigenvalues and eigenvectors we can 
now ask if our system of equations is hyperbolic. Con- 
sider first the case / = 0. In such a case one can show 
that the eigenvalues reduce to 



(62) 



that is, both eigenvalues are equal and the eigenvec- 
tors (61) do not form a complete set, so the system is 
only weakly hyperbolic. 

In the case when / < 0, one can use the fact that 
II 2 > * 2 > viz 2 (which will hold for spacelike slices) to 
show that the term inside the square root in equation (60) 
is always negative and the system is not hyperbolic. 

The / > case turns out to be more difficult to an- 
alyze. If < / < 1, then it is not difficult to prove that 
the term inside the square root is always positive. This 
means that we have two distinct real eigenvalues and a 
complete set of eigenvectors, and the system is therefore 
strongly hyperbolic. If / > 1, on the other hand, the 
term inside the square root can become negative for suf- 
ficiently large ^ 2 :— • ^ would then appear that 
in such a case we do not have a hyperbolic system. How- 
ever, the fact that hyperbolicity depends on the size of Vl/ 2 , 
indicates that this is only a coordinate problem. Indeed, 
if we reorient our spatial coordinates in a way such that 
only has components along the (t, x) directions, then 
'J 2 vanishes and the eigenvalues become real. The fact 
that for a different orientation of the spatial coordinates 
we can have complex eigenvalues is just an indication 
that for / > 1 the characteristic cones can tilt beyond 
the (t, x) coordinate plane. This is analogous to solving 
the hydrodynamic equations using a supersonic reference 
frame: the hydrodynamic equations become elliptic, but 
this is just a consequence of choosing a bad coordinate 
system. We then conclude that for / > 0, an orientation 
of the coordinates always exists such that our system of 
equations is strongly hyperbolic. 

Using the expression for the eigenvectors above, the 
condition for linear degeneracy, Eq. (45), takes the form 



d\± d\ ± 



(63) 



A straightforward calculation gives the following inde- 
pendent linear combinations of the previous conditions 



C+ + C*_ = = a 2 n [2a (1 + a) + aa'} 
( a 2 n 2 (1 + a + aa 2 * 2 ) + 3a 2 (IT 2 - * 2 ) - 3 ' 
\ (l + aa 2 n 2 ) 3 

C+ - C- = = -a 2 * 2 [2a (1 + a) + aa'} 
3a 2 II 2 (1 + a + aa 2 * 2 ) + a 2 (II 2 - * 2 ) - 1 
(1 + aa 2 n 2 ) 3 Jl + aa 2 * 2 



,(64) 



,(65) 



where a' := da/ da. For arbitrary a, II and the only 
way in which both these conditions can hold is if we take 



2a (1 + a) + aa' = 
which turns out to be equivalent to 

1 - / - af/2 = . 



(66) 



(67) 



This is precisely the condition (39) found in reference [6] . 
The main difference between the derivation above and the 
one of reference [6] is that in that reference the condition 
was derived using a hyperbolic formulation of the Ein- 
stein equations (the Bona-Masso formulation) , while the 
new derivation is based purely on analyzing the slicing 
condition and makes no reference to the Einstein equa- 
tions in any form, showing the generality of the condition. 

As already shown in reference [6], condition (39) can 
be trivially solved to find 



/ = 1 + k/a 2 



(68) 



with k > an arbitrary constant. For k — we recover 
harmonic slicing. On the other hand, if k ^ we see that 
for small a we have / <~ a~ 2 , so the results of Sec. Ill 
imply that the slicing will be strongly singularity avoid- 
ing. However, we can also see that in this case we are 
in the regime for which the lapse can easily become neg- 
ative. This means that the solution (68) has a serious 
drawback for any non-zero value of k, since it can allow 
the lapse to become negative as it collapses toward zero. 

In the next sections we will see how one can still ob- 
tain useful slicing conditions by looking for approximate 
solutions to condition (39). 



C. Zero order shock avoidance 

In the previous section we found that if we want to 
guarantee that no shocks will form, then we must choose 
the function f(a) in a way that is incompatible with hav- 
ing a lapse that does not become negative when it col- 
lapses toward zero (except in the specific case of harmonic 
slicing). Here we will relax our requirements and look for 
approximate solutions of condition (39). 

We start by assuming that the lapse is very close to 1, 
that is 



l + e 



(69) 



with e«l. Notice that the limit above applies to sit- 
uations that are close to flat space, but generally docs 
not apply to strong field regions (like the region close 
to a black hole) where the lapse can be expected to be 
considerably smaller than 1. However, in such regions 
considerations about singularity avoidance are probably 
more important. Our aim is to find slicing conditions that 
can avoid singularities in strong field regions, and at the 
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same time don't have a tendency to generate shocks in 
weak field regions. 

We can now expand / in terms of e as 

/ = ao + ai£ + 0(e 2 ) 

= aa + ai (a-l)+0(e 2 ) , 

and look for solutions to (39) to lowest order in e. 
Substituting (70) into (39) we find 



(70) 



a - ai /2 + O(e) = 



(71) 



This means that if we want condition (39) to be satisfied 
to zero order in e we must have 



oi = 2 (1 - a ) , 

which implies 

/ = a + 2(l-a Q )e + O(e 2 ) 

= (3d -2) + 2(l-d )a + O(e 2 ) . 



(72) 



(73) 



We must remember that (73) is just an expansion for 
small e. Any form of the function f(a) that has the same 
expansion to first order in e will also satisfy condition (39) 
to zero order. One family of such functions emerges if we 
ask for f(a) to have the following form 



Po 



(74) 



It is not difficult to show that for this to have an ex- 



pansion of the form (73) we must ask for 



p = a , gi = 2(d - l)/d , 



(75) 



which implies 



/ = 



a + 2 (a - 1) e 

<4 

(2 - d ) + 2 (ao - 1) a 



(76) 



Notice that if we take do = 1, we recover harmonic slic- 
ing. But there is one other case that is of special inter- 
est: For a = 2 the previous solution reduces to / = 2/a, 
which corresponds to a member of the 1+log family. The 
crucial observation here is that, as already mentioned in 
Sec. II A, this specific member of the 1+log family is pre- 
cisely the one that has been found empirically to be very 
robust in black hole simulations [13, 14, 24]. The fact 
that it is the only member of the 1+log family that sat- 
isfies condition (39) even approximately means that one 
should in fact expect it to be particularly well behaved. 



D. First order shock avoidance 

We can now go one order higher in e to obtain other 
interesting forms of /. Taking now 

/ = a + die + a 2 e 2 + C(e 3 ) 

= a + ai (a-l) + a 2 (a-l) 2 + 0(e 3 ) , (77) 



we find, after substituting into condition (39), that 

(1 - d - di/2) - (3di/2 + d 2 ) e + 0(e 2 ) = . (78) 
Condition (39) will be satisfied to first order if we take 



di = 2 (1 - d ) , 

d 2 = -3di/2 = -3 (1 - d ) 



(79) 
(80) 



So the expansion of / takes the form 

/ = do + 2(l-d )e-3(l-d )e 2 + 0(e 3 ) 
= (6d — 5) + 8 (1 — do) a 

-3(1 - d ) a 2 + C(e 3 ) . (81) 

Just as before, we can now look for rational functions 
that have the above expansion. One such possibility is 
to ask for / to have the form 



/ = 



Po 



(82) 



1 + die + q 2 e 2 

In order for / to have the expansion (81) we must take 

Po = d , (83) 
oi = 2 (d - 1) /d , (84) 
52 = (1 - do) [3d + (1 - d )] /a 2 , (85) 

which means that / takes the final form 



/ = 



a 2 , - 2d (1 - d ) e + (1 - d ) (4 - a ) e 2 



(4 - 3d ) + a (1 - d ) [(4 - d ) a - 8] 



(86) 



If we take d = 1 we again recover harmonic slicing. 
Another interesting case is obtained by asking for 



4 - 3d = d = 4/3 



(87) 



since in that case we will have / <~ oT 1 for small a, and 
as we have seen this implies good singularity avoidance. 
For such a choice the function / reduces to 



/ 



3a (3 — a) 



(88) 



For small a, this form of / behaves as a member of the 
1+log family. Moreover, it satisfies condition (39) to 
higher order than the usual choice / = 2/a. One could 
worry about the fact that for the choice (88) the func- 
tion / can become negative for a > 3. However, such 
a situation is unlikely to arise in practice since the ini- 
tial lapse is always taken to be at most 1 throughout the 
region of interest, and later evolution usually makes it 
even smaller. Because of these facts, the slicing condi- 
tion given by the choice (88) would seem to be a good 
candidate for a robust slicing condition when evolving 
systems with strong gravitational fields (including black 
holes). Whether this is true or not can only be settled 
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by numerical experiments where this form of / is tested 
against more traditional choices. 

Notice also that, whereas in the case when / = 2/ a 
the asymptotic gauge speed in regions where a ~ 1 is 
\/2 <~ 1.41, in the case (88) the gauge speed in those 
regions is only -y/4/3 ~ 1.15, which is much closer to 
the physical speed of light and might represent an extra 
advantage as gauge effects will not propagate much faster 
than physical effects. 

V. CONCLUSION 

I have considered the Bona-Masso hyperbolic family of 
slicing conditions and have studied under which circum- 
stances such hyperbolic slicings can avoid two different 
types of singularities: focusing singularities, defined as 
those for which the spatial volume elements vanish at 
a bounded rate, and gauge shocks, defined as coordinate 
singularities for which the lapse becomes discontinuous as 
a consequence of the crossing of the characteristic lines 
associated with the propagation of the gauge. 

In the case of focusing singularities, I have extended 
the analysis of Bona et. al. and shown that, depending 
on the form that the function f(a) defining the slicing 
takes in the limit of small a, one can have three differ- 
ent types of behavior: the lapse vanishes before the spa- 
tial volume elements do (collapse of the lapse and strong 
singularity avoidance), the lapse vanishes at the same 
time as the spatial volume elements and the singular- 
ity is reached after an infinite coordinate time (collapse 
of the lapse and marginal singularity avoidance), or the 
lapse vanishes at the same time as the volume elements 
but the singularity is still reached in a finite coordinate 
time (collapse of the lapse but no singularity avoidance). 
Harmonic slicing falls into the marginal singularity avoid- 
ing case, whereas the more commonly used 1+log family 
falls into the strong singularity avoiding case. 

For the case of gauge shocks I have re-derived, in a way 
that is independent of the Einstein equations, a condition 
on the function f(a) found previously that avoids them. 



This condition, unfortunately, is severely restrictive and 
implies that the lapse can easily become negative dur- 
ing evolutions. I have therefore studied different forms of 
the function /(a) that satisfy the condition only approx- 
imately. This study has shown that one specific mem- 
ber of the 1+log family that has previously been found 
empirically to be particularly robust, is in fact the only 
member of that family that satisfies the condition for 
shock avoidance even to lowest order. By asking for the 
shock avoidance condition to be satisfied to higher order, 
I have found a new form of the function f(a) that has 
the potential of being even more robust than the 1+log 
slicings for simulations of strongly gravitating systems. 

As a final comment, it is important to mention that el- 
liptic slicing conditions such as maximal slicing can eas- 
ily avoid both focusing singularities (since the volume 
elements are not allowed to change) and gauge shocks 
(since the speed of propagation is infinite), so they should 
in principle be more robust that the slicings considered 
here. Elliptic conditions, however, are considerably more 
computationally expensive. Not only that, but they are 
typically much more restrictive. For example, maximal 
slices might not always exist (they typically do not ex- 
ist in cosmological scenarios). Finally, elliptic equations 
require boundary conditions that might be difficult to 
impose in some situations, in particular when one has 
internal boundaries such as those associated with black 
hole excision. It is because of these reasons that one 
should consider alternatives to elliptic gauge conditions, 
such as those studied here. 
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